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The dynamics of quasicrystals is characterized by the existence of phason excitations in addition 
to the usual phonon modes. In order to investigate their interplay on an elementary level we resort 
to various one-dimensional model systems. The main observables are the static, the incoherent, 
and the coherent structure factor, which are extracted from molecular dynamics simulations. For 
the validation of the algorithms, results for the harmonic periodic chain are presented. We then 
study the Fibonacci chain with harmonic and anharmonic interaction potentials. In the dynamic 
Fibonacci chain neighboring atoms interact by double-well potentials allowing for phason flips. The 
difference between the structure factors of the dynamic and the harmonic Fibonacci chain lies in 
the temperature dependence of the phonon line width. If a bias is introduced in the well depth, 
dispersionless optic phonon bands split off. 
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I. INTRODUCTION 
A. Motivation 

Quasicrystals are long-range ordered materials lacking 
translational symmetry^. Their difTraction patterns ex- 
hibit a dense set of sharp Bragg reflections, that can 
be indexed by an integer linear combination of a finite 
number of basis vectors which is larger than the dimen- 
sion of space. As a consequence the atomic configuration 
of quasicrystals is described with reference to a higher- 
dimensional analog of a periodic lattice. Elementary dy- 
namic excitations within this 'hyperspace' description are 
phonons^ and phason^. 

Phasons involve rearrangements of the structure by 
atomic jumps over short distances, denoted 'phason 
flips'. They are connected with many physical proper- 
ties of quas icry stals as f or e xample elastic deformations^, 
dislocations^^, diffusion^!, and pha se tra nsformations^. 
Recently, indications for phason flip d 10 * 11 ! have been ob- 
served by in situ transmission electron microscopy. A 
coherent set of phason flips may form a static phason 
field, e.g. during a phase transformation or in the neigh- 
borhood of a dislocation^. 

By investigating the dynamics of quasicrystals one 
can find out the influence of the quasiperiodicity on the 
phonon spectrurrP and one may gain a deeper under- 
standing of phason flips 13 . Both points can be studied 
in x-ray or neutron diffraction experiments by measur- 
ing the response of the system in frequency (w) and 
momentum (q) space. Depending on the experimental 
setup, different functions can be obtained from the scat- 
tering experiments: (1) The static structure factor S(q) 
is the usual - not energy resolved - diffraction image, 
measured with cither x-rays or neutrons. It is used for 
the determination of the atomic structure. (2) The co- 
herent structure factorSS S(q,uj) is studied via coherent 



inelastic neutron scattering^ or alternatively via inelas- 
tic x-ray scattering^. It allows the determination of the 
phonon dispersion relations. The experiments on icosa- 
hedral quasicrystals show well defined acoustic phonon 
modes at small wave- vectors^ and dispersionless broad 
optic bands at larger wave vectors^. The cross-over be- 
tween the two regions is very sharp. (3) The incoherent 
structure factor^ Si(q, u>) can be measured in quasielas- 
tic neutron scattering. Neutrons are exclusively used 
here, due to the necessity of a high energy resolution. 
The technique also allows the investigation of phason 
flips. In a series of experiments Coddens et alfUMMEl] 
have found an anomalous (/-dependence of the quasielas- 
tic signal in icosahedral quasicrystals. They interpreted 
it as correlated simultaneous jumps of several atoms.^ 

Up to now various calculations of the coherent struc- 
ture factor of quasiperiodic model systems have been 
published, see Ref. |U Amon gst them are the perfect one- 
dimensional Fibonacci chain™! by static phason fields 
disordered Fibonacci chains^", and three-dimensional 
tilings^. In these studies the dynamical matrix is diago- 
nalized, which is a purely analytic method and yields the 
phonon dispersion relations only. The results are highly 
structured excitation spectra with a hierarchical system 
of gaps^l. The influence of anharmonicities, however, es- 
pecially the dynamics of phason flips has not been taken 
into account. 

This 'missing link' markes the starting point of our 
study. Here we present calculations of the structure fac- 
tors of special one-dimensional quasiperiodic model sys- 
tems by use of molecular dynamics (MD) simulations 
with either harmonic potentials or potentials that allow 
for phason flips. Although structure factors play such a 
central role in the dynamics of solids, not much seems 
to be known about their exact forms for one-dimensional 
chains. E ven for the simple harmonic chain only few ar- 
tides exispESEaSH. 
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FIG. 1: Double-well potential V(x) — x 4 — 2x 2 of the dynamic 
Fibonacci chain. The equilibrium distances are S and L. 

B. Model systems 

As a simple one-dimensional model for a quasiperiodic 
system we consider the Fibonacci chain. It consists of 
particles arranged at two different distances: large ones 
(L) and small ones (S). The length ratio L/S equals 
the number of the golden mean r = 1). The 

sequence of the distances is created recursively by the 
mapping {L, S} i— > {LS, L} with starting condition L. 
For example, after four iterations the resulting sequence 
is LSLLSLSL. 

We want to study chains consisting of identical parti- 
cles with nearest-neighbor interactions. The Hamiltonian 
has the form 

N 2 

« = Ei + ^r%i-4 (1) 

3=1 

where Xj and pj are position and momentum of the jth 
particle. The dynamic Fibonacci chain (DFC) is defined 
by the choices aj = an = r 3 for the equilibrium distances 
and V(x) — x 4 — 2x 2 for the interaction potential, re- 
spectively. The latter forms a double-well potential with 
minima at ±1 and a potential hill of height AE = 1 
as shown in Fig. [I] Because neighboring particles sit in 
either of the potential minima, two nearest-neighbor dis- 
tances L = a + 1 and S = a — 1 are possible. They 
fulfil the constraint L/S = t of the Fibonacci chain. 

The DFC shows two types of elementary excitations: 
Phonon vibrations in the local minima and phason flips 
that interchange the particle distances L and S. At low 
temperatures only phonons are excited, phason flips have 
to be activated thermally. With its neighbors at rest the 
activation energy of a particle for a phason flip is 2AE. 
This value is a result of the perfect superposition of the 
potential hills of both neigbors. The value is lowered 
when the neighbors assist by stepping simultaneously to 
the inside or outside during the phason flip thus creating 
a non-perfect superposition of the potential hills. Since 
the particle distances L and S are energetically degen- 
erate, the total equilibrium potential energy is invariant 
under a phason flip. 

The occurrence of phason flips makes nonlinearity an 
intrinsic feature of the DFC and an analytical treatment 



of the dynamics impossible. To understand the influence 
of the nonlinearity, we study four model systems with 
increasing complexity concerning their dynamical behav- 
ior: 

• Harmonic periodic chain (HPC): 
Vhpc(x) — 4a; 2 and aj = a = 2\/5. 

• Harmonic Fibonacci chain (HFC): 

Vhfc( x ) — 4x 2 and aj = L or S according to the 
Fibonacci sequence. 

• Dynamic Fibonacci chain (DFC): 
Vdfc( x ) = x A — 2x 2 and aj = t 3 . 

• Asymmetric Fibonacci chain (AFC): 

Vafc(x) = l/ D FC + x(^ 2 -l) 2 (« + a; 2 /2-l/2) and 
aj = t 3 with x € [0, 1] and e = ±1. 

The potentials of the HPC, HFC, and DFC are chosen to 
be identical in the harmonic approximation around the 
equilibrium separation. The average particle distance a is 
the same for all four systems. In the case of the Fibonacci 
chain the occurrence probabilities for L and S are given 
by t -1 and t~ 2 , hence a = 2y5- 

For the solution of the equations of motions we use 
a special molecular dynamics (MD) code. The code is 
introduced in Sec. |TT] together with a short theoretical 
background. The simplest system is, of course, the HPC. 
Exact solutions for the equation of motion exist as a su- 
perposition of plane waves. We study the dynamics of the 
HPC in Sec. |III| as a reference system. The HFC con- 
sists of particles arranged on the Fibonacci chain with 
distances L and S interacting with the same harmonic 
potentials as the HPC, see Sec. [IV] The DFC will then 
be studied in Sec. |Vj In the case of the AFC the particles 
in the two potential minima have different eigenfrcquen- 
cies. The parameters \ and e determine the degree of 
asymmetry. For more details we refer to Se c. |VI| We 
finish with a discussion and conclusion in Sec. I VIII 

II. STRUCTURE FACTORS FROM 
MOLECULAR DYNAMICS 

A. Definition of the structure factors 

We write the particle number density of the chain with 
N particles as a sum of delta functions positioned along 
the particle trajectories xi(t), n(x,t) = Y^iLi S(x — xi(t)). 
The time dependent density-density correlation function 
and the density-density autocorrelation function are de- 
fined as 

G(x,t) = ^ J (n(x',t)n(x + x',0)) dx' 

= iE^-^W + ^CO))}, (2a) 

3,1 

G*(x,t) = ~^<£(z-x;(*) + ^(fj))>, (2b) 
I 
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where the brackets denote the thermal average P3 The 
coherent and incoherent structure factor are the space- 
time Fourier transforms^!, 



and 



S(q,u>) = ^-^ J e- iuJt J2{ e ~ lqX3it)eiqXl{0) ) dt > 

3,1 



Si(q,cj) = — "-r f e- luJt J2( e ~' lqXl{t)elqM0) ) dt - 



(3a) 



(3b) 



Both functions are symmetric about q — and ui = 0. 
The static structure factor is the integral of the coher- 
ent structure factor, S(q) = J S(q, to) du>, i.e. the Fourier 
transform of G(x, 0), 



3,1 



e -iqxj(0) e iqxi(0) 



(4) 



B. Molecular dynamics simulations 

For further calculations the particle trajectories are re- 
quired as solutions of the equations of motion. Since in 
the case of the anharmonic chains only numerical solu- 
tions exist, we use a simple MD code. Initially the par- 
ticles are placed on the equilibrium positions of a finite 
chain of length L with periodic boundary conditions. The 
velocities are initialized according to a Gaussian distribu- 
tion. Its width determines the total energy and thus the 
temperature of the system. The equations of motion are 
integrated by a Verlet algorithm running for a simulation 
time T s j m . After starting the simulation, the dynamics is 
not controlled by a thermostat or in any other way. 

For the direct numerical calculation of the Eqs. ^ we 
must compute a fourfold sum: two sums over the par- 
ticle number N and two over the time T s j m , one sum 
for the Fourier transform and one for the time average. 
Note, that by assuming ergodicity the thermal average 
( ) can be replaced by a time average j^— J dt and ad- 
ditionally by an average over several independent MD 
runs. For the sake of clarity the averaging over the MD 
runs is suppressed in the following notation. We intro- 
duce a more compact notation by defining the functions 
fi(q,t) — e lqxi ^\ Let us assume tentatively that these 
functions are periodic in time with period T B i m and in 
space with period L. Then the Eqs. (|3| and Q are 
greatly simplified to 



S{q,w) = 



1 



J2fl(q,t)dt 



, (5a) 
(5b) 



S(q) = 



NT si 



dt. 



(6) 



The equations differ in the order of the absolute square 
and the particle sum. Since only two sums are left, an 
efficient numerical computation of the structure factors 
is possible. Furthermore a fast Fourier transform is used 
for the time integrals in Eqs. ([5|. 

It is left to discuss the periodicity conditions. The spa- 
tial periodicity follows from the periodic boundaries used 
in the simulation. Therefore the chain acts like a ring and 
the excitations can go round during the simulation. To 
avoid such a behavior that would lead to unwanted cor- 
relations, we limit the maximum simulation time by the 
quotient of the length of the chain and the sound velocity 
c s to T max = L/c s . For all the model systems HPC, HFC, 
and DFC the sound velocity is the same, c s = a\/8, and 
T max = N/y/8. We use T slm = T max . 

The functions fi(q,t) are in general not periodic in 
time. There is no reason, why the particles should be 
at the same positions at the end of the simulation as 
at the beginning. To avoid this problem, we multiply 
fi(q,t) with a window function w(t) to enforce an arti- 
ficial periodicity. The function w(t) has to decrease fast 
enough - both in direct as in Fourier space - towards the 
boundaries of its domains. We use a normalized broad 
Gaussian function. Its width is chosen as large as pos- 
sible with the constraint that the Gaussian has decayed 
to a small enough value at the interval boundaries. The 
effect of the Gaussian is a smoothing of the structure fac- 
tors by convolution with a narrow Gaussian. The exact 
value of the width has no influence on the results. 



III. HARMONIC PERIODIC CHAIN AS 
REFERENCE SYSTEM 

A. Analytic calculations 

The harmonic periodic chain (HPC) is used as a ref- 
erence system to test our algorithms since its equa- 
tions of motions can be solved analytically. If we put 
xi(t) = ui(t) + la, then the ui(t) are expressed by a linear 
combination of normal modes. The wave vector q and 
the frequency u> are related according to the dispersion 
relation u(q) = 2tt;o | sin(ga/2) | . Here, wq is the eigen- 
frequency of a single particle. In the case of the model 
systems HPC, HFC, and DFC we have w = 

For the HPC the thermal averages in the structure 
factors, Eqs. Q can be calculated to be 

oo 

S(q,cj) = - J 4 e-^exp(-ig 2 C rf(t))dt 1 



SiM = 2^f-^2 je- iut fi(q,t)dt , (5b) = ±_ j 6 xp (-igV 2 (t)) dt. 



(7a) 
(7b) 
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where we used from the literature^ 



*?(t) 



(Mi) 
k B T 



«o(0)] 2 ) 



2u t 



J2i(s)(2u>Qt — s) ds 



(8) 



J n (s) is the Bessel function of the first kind of order 
n. The only external parameter in these functions is 
the temperature T. The particle sum and the Fourier 
transform have to be evaluated numerically. Due to the 
translation invariance of the HPC the double sum of the 
Eqs. ^ is reduced to a single sum. The static structure 
factor of the HPC can also be calculated analytically 26 



S(q) 



sinh(q 2 cr 2 /2) 



cosh(<7 2 cr 2 /2) — cos(ga) ' 



(9) 



where a 2 = k B T/tu 2 . 

Let us take a closer look at the incoherent struc- 
ture factor Si(q,u>). In the limit of small T the term 
cxp(— \q 2 a 2 {t)) decays slowly with t and we substi- 
tute ao(t) with its approximation for large t: ao(t) = 
\t\kBT/ojQ for \t\ — > oo. This leads to a Lorentzian peak 



1 



7r r 2 



r 



r 



q 2 k B T 
2cu 



(10) 



In the limit of large T the term ao(t) is approximated for 



small t: <Jo(t) = 
Gaussian peak 

Si(q,w) 



t 2 k B T for \t\ < 



Hence there is a 



T 



'2n 



exp 



2i- 



r), 7 = Vq 2 k B T. (11) 



The transition temperature between these two limiting 
cases is k B T = 4w 2 /<? 2 . 

The Fourier transform of the Eqs. Q yields the corre- 
lation functions 



G(x,t) 



G a (x,t) 



l y> l 



l=-c 



exp 



[x + la)' 



/2ncr (t) 



exp 



2a 2 {t) 



(12a) 
(12b) 



The function G(x,t) is shown in Fig. [2j It consists of a 
sum of Gaussians centered at the equilibrium positions 
of the particles. The width of the Gaussians increases 
with temperature, as well as with time and in space: 
lim a . ! t_> 00 G(x, t) — p and p — I /a. This means that 
there is no longe-range order. Indeed, for the particle 
number density we have (n(x)) = p. which is uniform as 
in liquids^. The autocorrelation function G a (x, t) corre- 
sponds to the center peak at x = 0. 



B. Simulation results 

In the case of the incoherent structure factor the vari- 
ables T and q only appear as combination q 2 k B T in 




FIG. 2: Density-density correlation function G(x, t) of the 
HPC for ksT — 0.5. The Gaussians are centered at integer 
multiples of x — a w 4.47. 
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FIG. 3: Incoherent structure factor Si{q,uj) of the HPC for 
q — n /a, and different temperatures. The symbols mark the 
data from MD simulations with iV = 1000 p arti cles, the lines 
the result from the analytical fomula, Eq. ( |7b[ ). The peaks 
and edges are at integer multiples of u> = 2u>o = 2v / 8- 



Eqs. ^ and Therefore it suffices to examine Si at a 
fixed wave vector for different temperatures. We choose 
q = tt / a arbitrarily. The results from MD simulation 
and the numerical integration of the analytical formula, 
Eq. (7b I are shown in Fig. [3] for temperatures ranging 
from 0.01 to 100.0. There is a maximum at lu = 0, called 
the quasielastic peak. In the low temperature regime the 
maximum has a Lorentzian shape. Furthermore a one- 
phonon peak, a two-phonon edge, a three-phonon edge, 
etc. are found at 2ljq w 5.7, 4cj ~ H-3, ~ 17. 0.^ 
Note that the multi-phonon contributions rapidly decay 
for larger oj (logarithmic scale). At higher temperatures 
the curve smoothes and approaches a Gaussian profile. 

The coherent structure factor S(q, lu) is shown in Fig. 
[4] For a low temperature of k B T = 0.01 a one-phonon 
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branch, a two-phonon branch, and very weakly a three- 
phonon branch are observed. The one-phonon branch 
has a Lorentzian line shape and follows the phonon- 
dispersion relation. At higher ^-values the branch broad- 
ens with a width proportional to ksTq 2 . This is similar 
to the temperature behavior for the incoherent structure 
factor. The multi-phonon branches follow the modified 
relations to(q) = 2nuJo\ sin(qa/2n)| with n = 2, 3. 

In Fig. Qb) the comparison of the MD simulation (left 
side) and the analytical formula Eq. (7a I (right side) is 
shown. The temperature in this figure is kgT = 0.1, 
which is higher than in Fig. |3Ja). As a consequence 
the one-phonon branch is broader. For both methods 
of calculating S(q,u>) a high accuracy over 12 orders of 
magnitude is possible. The accuracy is only limited by 
the internal floating point precision. 

The static structure factor from Eq. ^ is compared 
to the results from the MD simulation in Fig. [5j S(q) 
consists of a sequence of Lorentzian peaks at the recipro- 
cal lattice points. The increasing width for larger wave 
vectors shows again that no long-range order is present 
in the one-dimensional model system. 

The MD simulations and analytical formulas show a 
perfect agreement. Therefore we conclude that MD sim- 
ulations are a well-suited numerical tool for calculating 
the structure factors of the one-dimensional model sys- 
tems. Although we integrate the equations of motions 
with a good precision only on a short time scale using the 
simple Verlet-algorithm, the statistics extracted from the 
trajectories are correct. This confirms our approach and 
encourages us to proceed studying the phason dynamics 
of Fibonacci chains in the next section. 



IV. HARMONIC FIBONACCI CHAIN: 
INFLUENCE OF THE QUASIPERIODICITY 

By changing the interparticle equilibrium distances of 
the HPC to those of a Fibonacci sequence with separa- 
tions L and S we obtain the harmonic Fibonacci chain 
(HFC). The interaction potential is left unchanged. Since 
the incoherent structure factor is a function of the single 
particle motion only, it does not depend on the equilib- 
rium distances of the particles but only the interaction 
potential. Hence, the incoherent structure factor of the 
HPC and of the HFC are identical. For the coherent 
structure factor the interparticle distances become im- 
portant. Instead of Eq. (|7a|) we now have 



S(q,w) 



1 

27 



Jt E< 



-%qx\ 



exp i 



„o _ 



Here x° = Y,\=i a j for 1 > 0; x< { = 



q 2 af(t)/2) dt. 

(13) 

for I < and 



=/ 3 



denote the equilibrium positions of the particles. 
For UbT = this gives the Fourier transform of the static 
Fibonacci chain 



S(q,uj) = 5(cj)^V 



(14) 



which is well knowrH It consists of a dense set of delta 
peaks with varying intensity, positioned at the reciprocal 
lattice points 



q = q Q {h + Th!) , 



(15) 



with g = ^ « 1.40. 

As shown in Fig. [6] for the temperature ksT = 0.02, 
the coherent structure factor S(q, u>) of the HFC consists 
of many different branches all following the one-phonon 
dispersion relation. Relative to each other the branches 
are displaced in the q direction. They start at the recip- 
rocal lattice points with the intensity of the respective 
delta peak. Two-phonon branches are also found. The 
broadening of the branches proportional to fc^T^ 2 has 
already been discussed for the HPC. 

From these findings and Fig. [5] one can guess the 
form of the static structure factor S(q) for the HFC: 
Lorentzians are positioned at the reciprocal lattice points 
with varying intensity. In Fig. [7] S(q) is shown for 
q £ [0,37r]. For small q a large number of Lorentzians 
occur. With increasing q-value their width increases and 
the stronger ones hide the weaker ones. Comparing the 
position of the strong peaks and ignoring the change of 
their widths and heights, there is a self-similarity in S(q). 
The deflation factor is Aqi/Aq 2 — r 3 as indicated in 
Fig. [7] 

It is interesting to note that there are regions with very 
few peaks. They are positioned around q — tt, q = 2ir, 
etc. . The same regions are also special for the coherent 
structure factor. As seen in Fig. [6] all the weak one- 
phonon branches vanish towards q — n. Only the strong 
one-phonon branch starting from the Bragg peak at q = 
rqo ss 2.27 remains. 



V. DYNAMIC FIBONACCI CHAIN: 
OCCURRENCE OF PHASON FLIPS 

A. Phason flips 

Let us now proceed to the anharmonic chains with pha- 
son flips by investigating the dynamic Fibonacci chain 
(DFC). It is built from identical particles that interact 



with a symmetric double-well potential Vdfc 



2x 2 



First the notion of a phason flip has to be specified. To 
do so we identify the position of the changes from L to 
S and from S to L of the interparticle distances along 
the particle trajectories. This is done in Fig. [8] Often 
a L — > S change and a S — > L change lie next to each 
other (particle distance 1) and the sequences LS and SL 
are interchanged. But there are also many cases where 
the positions of the two changes are separated by 0, 2, 
3, or even more particle distances as marked by lines in 
Fig. [8] Sometimes it is not possible to find a partner lo- 
cally. Only in the long time average every L — > S change 
will eventually cancel with a S — > L change. 

In the literature on the Fibonacci chain a phason flip 
is understood as the exchange of a L and a neighboring 



6 




<r 4 3 




FIG. 4: (Color online) Coherent structure factor S(q, w) of the HPC with N — 6500 particles from MD simulation. The 
temperatures are for ksT = 0.01 (a) and for ksT = 0.1 (b). One-, two- and three-phonon branches are observed. They start 
at the reciprocal lattice points 2nn/a. In (b) the output from the MD simulation (left side) is compared to the output of the 
analytical formula Eq. (7a I (right side). 




FIG. 5: Static structure factor S(q) of the HPC for k B T = 
0.02. The symbols mark the data from a MD simulation with 
N = 1000 particles, the line the result from the analytical 
fomula, Eq. j9j. 



S particle distance. As we learned above there are also 
other types of exchanges in the DFC. In the following we 
denote by phason flip every pair of flip partners as those 
connected by lines in Fig. [8] Note that the times and 
positions of the phason flips are not well defined. Only 
their number can be estimated by counting the changes 
in the particle distances as we will do now. 

The temperature dependence of the average phason flip 
frequency ojfa p is shown in Fig. [9] Phason flips start to 
appear at about k^T = 0.1. At low kgT the average pha- 
son flip frequency increases rapidly by thermal excitation 
and Wflip < wj. At higher temperatures k B T > 0.4 the 
average phason flip frequency slowly saturates. In this 
region the internal energy is comparable to the potential 
hill. 



B. Results for the structure factors 

For the anharmonic chains no analytic results are avail- 
able, in particular not at elevated temperatures when 
phason flips occur. The incoherent structure factors of 
the DFC and of the HFC /HPC differ remarkably and are 
shown in Fig. 10 for q = n/a. The comparison leads to 



the following conclusions: 

1. At a fixed temperature, there are w-ranges where 
the curve for the DFC lies below the curve for 
the HFC/HPC and vice versa. Since we have 
/ Si(q,uj) duj = 1 from Eq. (3b I, the integral area 
between the two curves has to vanish. 

2. At very low k B T and u < 2u)q ~ 5.7 the curves 
of the DFC and the HFC/HPC cannot be distin- 
guished in logarithmic scale except for two small 
bumps. They are a consequence of the anharmonic- 
ity of the interaction potential of the DFC and not 
related to the phason flips. At larger ui values the 
multi-phonon edges have different heights. 

3. Above k B T — 0.1 the one-phonon peak and the 
multi-phonon edges in the curves for the DFC 
broaden and weaken considerably faster than in the 
curves for the HFC/HPC. They finally disappear at 
k B T= 1.0. 

4. No additional peaks or edges occur at any temper- 
ature. 

Further MD runs show that different g-values change 
the temperature dependence, but the general features re- 
main unchanged: The phonon peaks and edges broaden 
and weaken much faster with increasing temperature for 
the DFC than for the HFC/HPC. 

Similar conclusions follow for the coherent structure 
factor of the DFC. S(q, u>) for the DFC looks qualitatively 
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FIG. 6: (Color online) Coherent structure factor S(q,uj) of the HFC with TV = 13000 particles for fcsT = 0.02. It consists of a 
dense set of phonon branches starting from the reciprocal lattice points. 
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FIG. 7: Static structure factor S(q) of the HFC from a MD simulation with 2000 particles at the temperature ksT = 0.02. 



similar to S(q,uj) for the HFC except that the branches 
broaden more quickly with temperature. To compare the 
broadening let us look at plane cuts through S(q, u>) for 
lo = 1.0 fixed. Three cuts for the temperature ksT = 
0.05, 0.2, and 0.3 are shown in Fig. [TT] A one-phonon 
branch is located inside the cut resulting in a peak with 
approximate Lorentzian line shape as indicated by the 
fits in the figure. Only for lower temperatures the line 
shape deviates from a Lorentzian, which is seen at the 
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is shown as 
121 for both, the HFC 



base of the peak for fc^T = 0.05. 

The width of the Lorentzian in Fig 
a function of temperature in Fig 
and the DFC. In the case of the HFC the width increases 
linearly with temperature, as discussed in Sec. [IVJ There 
is a deviation from linearity at low temperatures. This 
is an artefact from the method of calculation from MD 
simulation data. As explained in Sec. |IIB| the structure 
factor is convoluted with a Gaussian due to the finite 
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time 



FIG. 8: Snapshot of the particle trajectories of the DFC at 
the temperature ksT = 0.6. Changes in the particle distances 
from L to S and S to L are marked with a cross (x) and a 
plus (+). A L — > 5* change and a S — > L change combined 
form a phason flip. 
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FIG. 10: Incoherent structure factor Si(q,u>) of the DFC 
(solid) and the HFC/HPC (dashed) for q = 7r/a at different 
temperatures. The data was calculated using MD simulations 
with N = 1000 particles. 




FIG. 9: Average flip frequency as a function of the temper- 
ature fcgT. In the temperature range of the figure: u)mp < 
luo ~ 2.83. 
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FIG. 11: Cuts through S(q,u) of the DFC for a fixed U! = 1.0 
including a one-phonon peak. The solid curves are fits with 
Lorentzians. 



simulation time, 
observed offset. 



The (narrow) Gaussian generates the 
The convolution with the Gaussian is 

at 



also responsible for the shape of the cuts in Fig. 11 



low temperature, deviating from the Lorentzian shape. 

There is one aspect of the DFC that has been ignored 
up to now. Due to the phason flips the original per- 
fect quasiperiodic long-range order is slowly decaying. 
With progressing simulation time the chain becomes ran- 
domized which, however, has no effect on the incoherent 
structure factor. To test the influence of the randomiza- 
tion on the coherent structure factor, a MD simulation 
was started with a totally randomized Fibonacci chain. 
The interaction potentials and the occurrence ratio of L 
and S were not adapted. The result of the simulation is 
shown in Fig. 13 Most of the details in S(q, u>) are lost by 



the randomization process and the branches are strongly 



broadened, although the simulation has been carried out 
at a low temperature of UbT = 0.02. 

We summarize the results of this section: By the in- 
troduction of the anharmonic double-well potential of 
the DFC the phonon peaks, edges, and branches are 
strongly broadened and weakened with increasing tem- 
perature. There are two effects responsible for the broad- 
ening: (1) The anharmonicity and single phason flips 
affect the propagation of phonons. This is seen in the 
incoherent structure factor. (2) The destruction of the 
long-range order by a large number of phason flips. This 
is the main effect of broadening for the coherent structure 
factor. The phason flips appear more or less uniformly 
distributed over the simulation time and the particles. 
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k R T 

FIG. 12: Width of the Lorentzian peak of Fig. [TT] as a function 
of temperature for the HFC (inset) and the DFC. 




FIG. 13: (Color online) Coherent structure factor S(q, to) of 
a randomized Fibonacci chain with 2000 particles at the tem- 
perature ksT — 0.02. The same interaction potentials as in 
Fig. [6] have been used. 
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FIG. 14: Interaction potential Vafc(^) of the AFC for differ- 
ent values of the parameter \ an< A e = 1. Vdfc is shown with 
dashed lines. 



In harmonic approximation a particle will feel different 
eigcnfrequencies, depending on the nearest-neighbor con- 
figuration SS, SL, LS, or LL. In the case of e = 1 
it is uj\Aj > ^ls = <^sl > <^ss and for e = — 1 it is 
wll < ^ls = w sl < ^ss- The sign change of e corre- 
sponds to a mirror operation of Vafc about the y-axis. 

The coherent structure factor of the AFC with x = 0.1, 
0.3, 0.5 and e 



±1 is shown in Fig. 15 Band gaps of 
different widths appear and broaden with increasing val- 
ues for x- They are positioned at the frequencies where 
one-phonon branches intersect each other. Their posi- 
tions and widths are different for e = 1 and e = — 1. 
For e = l three large gaps and several smaller gaps ap- 
pear, whereas for e = — 1 only one very large gap, one 
medium gap and several small gaps appear. The band 
gaps are a consequence of the competing eigcnfrequen- 
cies due to the asymmetric potential in the same way as 
band gaps appear in periodic systems with several parti- 
cles and eigenfrequencies per unit cell. 



VI. ASYMMETRIC FIBONACCI CHAIN: 
COMPETING EIGENFREQUENCIES 



B. Density of states 



A. Band gaps 

In the next step we modify the double well potential of 
the DFC. Like all the previous model systems the asym- 
metric Fibonacci chain (AFC) is built of identical par- 
ticles, but they interact with the more complicated po- 
tential Vafc = Vdfc + AV, see Fig. 14 The additional 
term is given by 



AV(x) = X (x 2 - l) 2 (ex + x 2 /2 - 1/2) 



(16) 



with x € [0,1] and e = ±1. This term has been cho- 
sen in such a way that the positions of the minima at 
x — ±1 are left invariant, but the curvatures of the poten- 
tial at those points is changed to V£' FC (±1) = 8(1 ± ex)- 



The gaps are seen more clearly in the density of states 
(DOS) D(oS). The DOS per particle is calculated from 
the velocity autocorrelation function by a Fourier trans- 
form 



D(u) = 



ttN 



Y.iMtumdt, co>o. (17) 



It is obtained from MD simulation data similar to the 
incoherent structure factor: fi(q, i) has to be substited by 
vi (t) in Eq. ( |5b[ ) . By interchanging the Fourier transform 
and two time derivatives we can alternatively write 

D(LU) = £ fa(*M0)) dt - 




FIG. 15: (Color online) Coherent structure factor S(q, uS) of the AFC for different values of \ an< l e from MD simulations with 
6500 particles at ksT = 0.01. Band gaps are observed. The band gaps also appear in the DOS D(lj). 
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After a taylor expansion of the exponentials in Eq. ( 3b ) 



a connection to the incoherent structure factor is found. 



D(uj) = 2uj 2 lim 



Si(q,u) 



(19) 



For comparison we present the DOS of a harmonic 
chain (HPC or HFC) 



D(u) 



7T y^g - W 2 



for < lu < 2u! 



(20) 



and elsewhere. The total number of states is normalized 
to J °°L>(w)dw = k B T. 

It can be checked, that the DOS of the harmonic chain 
fits well to the DOS of the AFC for e = 0.0, except 
small bumps that originate from the anharmonicity of the 
potentials. In Fig. [l| the DOS of the AFC for different 
values of x and e are drawn. The band gaps appear at the 
same frequencies and the same widths as in the coherent 
structure factor. 



VII. DISCUSSION AND CONCLUSION 

At the end we would like to make some general re- 
marks concerning phason flips and phason modes: In the 
context of a hydrodynamic theory, phason flips can be as- 
sociated with phason modes. It was noted quite early^, 
that phason modes are diffusive, in contrary to the prop- 
agating phonons. This means that phason flips are only 
weakly coherent in space and in time, which of course we 
have also observed here. As a result their influence on the 
structure factor is small making it difficult though still 
interesting to study them by scattering experiments. An- 
other point concerns the connection of phason flips and 



quasiperiodicity. There is no reason why phason flips 
should only occur in quasicrystals. Since interaction po- 
tentials are not sensitive on the long-range order, phason 
flips in the form of atomic jumps can also occur in pe- 
riodic complex intermctallic p hases , which is supported 
by recent experimental results^^H. In the case of our 
model systems, it is equally possible to compare simula- 
tions of a periodic LSLSLS . . . chain with harmonic and 
double-well potentials respectively. 

In conclusion, we have investigated the dynamics of 
phonons and phason flips in one-dimensional model sys- 
tems with molecular dynamics simulations. An efficient 
algorithm made it possible to calculate the structure fac- 
tors with high precision and in great detail. As a result 
multi-phonon contributions, - although weak in compar- 
ison to the one-phonon peaks and branches - have been 
identified. By introducing phasons in the model systems 
we were able to study their influence on the structure 
factors, which is mainly a broadening of the character- 
istic peaks, edges, and branches with temperature. The 
broadening can be further split into a broadening due to 
the disorder as a result of collective phason flips, i.e. a 
static phason field, and a broadening due to the anhar- 
monicity of the interaction potentials and single phason 
flips. The work presented here is a first step. Further 
studies on two-dimensional and three-dimensional model 
systems with phason flips are under way. 
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